A walkthrough of Ernest (2005)’s original analytical approach, from close reading of the paper.
Ernest drew data from the Andrews LTER, the Sevilleta, Niwot Ridge, and Portal.
The data available online do not quite match the descriptive statistics reported in Ernest (2005).
replicate-becsDownload raw data. By default data will be stored in subdirectories of replicate-becs/data/paper/raw/ for each site.
download_raw_paper_data()
Process raw data into the appropriate format. This is a data table with a record for each individual and columns for species and weight in grams. By default these tables will be stored in subdirectores of replicate-becs/data/paper/processed.
process_raw_data()
Loading in data version 1.106.0
[1] TRUE
Load data tables for each community. There should be 9 communities.
communities <- load_paper_data()
length(communities)
[1] 9
Each community should be a data table with columns for species and size for each individual, for example:
names(communities)
[1] "andrews" "niwot" "portal" "sev-5pgrass" "sev-5plarrea"
[6] "sev-goatdraw" "sev-rsgrass" "sev-rslarrea" "sev-two22"
head(communities[[1]])
Ernest 2005 Fig 1
replicate-becsFor every individual, calculate metabolic rate and assign to a size class.
communities_energy <- lapply(communities, FUN = make_community_table, ln_units = 0.2)
head(communities_energy[[1]])
For each community, sum total energy use for each size class, and convert to the proportion of total energy use for that community.
bseds <- lapply(communities_energy, FUN = make_bsed)
head(bseds[[1]])
replicate-becsCalculate mean mass of each species in each community.
bsds <- lapply(communities, FUN = make_bsd)
head(bsds[[1]])
replicate-becsenergetic_dom <- lapply(communities_energy, FUN = energetic_dominance)
head(energetic_dom[[1]])
energetic_dom_prop <- lapply(communities_energy, FUN = energetic_dominance, mode_cutoff = 'prop')
replicate-becsbsed_uniform_bootstraps <- lapply(communities, FUN = community_bootstrap, bootstrap_function = 'bootstrap_unif_bsed_doi', nbootstraps = 10000)
See issue #4 on github.
replicate-becscommunity_combination_indices = utils::combn(x = c(1:9), m = 2, simplify = TRUE) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(community_a = V1, community_b = V2)
combine_communities = function(indices, communities) {
community_combination = list(community_a = communities[[indices[1]]], community_b = communities[[indices[2]]], community_names = c(names(communities)[[indices[1]]], names(communities)[[indices[2]]]))
return(community_combination)
}
community_combinations = apply(community_combination_indices, MARGIN = 1, FUN = combine_communities, communities = communities)
bsed_crosscomm_bootstraps = lapply(community_combinations, FUN = community_bootstrap,
bootstrap_function = 'bootstrap_crosscomm_bseds', nbootstraps = 10000)
See histogram of p values for comparisons to see if commuities’ BSEDs are the same or different.
Ernest (2005) refers to the communities with the site column above. To compare the communities above to the communities in the resurrected data set, we can try to match them based on the BSD and BSED plots (above) and species richness.
ernest_richness = read.csv(file.path("ernest-2005-files/ernest_richness.csv"), stringsAsFactors = F)
# Guesses based on body size plots
ernest_richness$community_name <- c('andrews', 'niwot', 'portal', 'sev-5pgrass',
'sev-rsgrass', 'sev-two22', 'sev-goatdraw',
'sev-5plarrea', 'sev-rslarrea')
bsds_richness = data.frame(community_name = names(bsds), stringsAsFactors = F)
bsds_richness$new_richness = NA
for(i in 1:nrow(bsds_richness)){
bsds_richness$new_richness[i] = nrow(bsds[[i]])
}
# See if the richness values match
ernest_richness = ernest_richness %>%
dplyr::left_join(bsds_richness, by = 'community_name') %>%
dplyr::mutate(richness_match = (richness == new_richness))
print(ernest_richness)
We can get 2 (of 7) communities at the Sevilleta to match up, and Niwot and Andrews. The other pairings are as close as possible.
Moving forward, we can compare the results of the Kolmogorov-Smirnov tests based on the values in the appendix and these name pairings.
ernest_key = ernest_richness %>%
dplyr::select(site, community_name)
write.csv(ernest_key, file = file.path(storagepath, 'ernest-2005-files/ernest_key.csv'), row.names = F)
ernest_bsds_uniform_results = read.csv(file.path(storagepath, 'ernest-2005-files/ernest_appendixA.csv'), stringsAsFactors = F) %>%
dplyr::left_join(ernest_key, by = 'site')
print(ernest_bsds_uniform_results)
replicate-becs:From Zar (1999) Biostatistical Analysis.
Tables of critical values were entered by hand from the appendix to Zar (1999).
The \(\delta\) corrected KS test does not correspond to the results from Ernest when the species mean body size values are on an untransformed scale.
Using the natural log of the species mean body size value, however…:
With mean mass logged, all the results replicate qualitatively (i.e. not significantly different from uniform) and Niwot, for which the currently-available data most closely matches that reported in Ernest (2005), replicates almost exactly numerically.
Ernest (2005) used a two-sample Kolmogorov-Smirnov test to compare every possible combination of community-level BSDs.
replicate-becs# use same community combinations as before
bsd_crosscomm_ks = lapply(community_combinations, FUN = ks_bsd,
ln_mass_vals = F)